/*
** Modified by Neal May/June 2011
** Remodified by Neal Sept 2014 for Letters piece
** Stripped out lots of things to just compare Logit and clogit for standard normal x
** and fixed earlier version to compute RMS error and bias
** This program is a modified version of
** 
** Bias in Conditional and Unconditional Maximum Likelihood Estimation
**  Ethan Katz (ekatz@fas.harvard.edu)
**  last modified 02/06/01
** of Tom Coupe's
** date: 02 02 2005
** this modification solves the error in Katz's original program as it does allow to include fixed effects
** Tom Coup?  tcoupe@eerc.kiev.ua
** NB - 2014 The modification is now almost 100% so I am no longer calling this a modification 
** Modified again for SPM 2015 presentation June 2015
** Note: Particular parameter values reflect whatever was the last run of the simulation
** To reproduce Table 1 make sure to reset the parameter values to those in the table
** Also note there are many things computed that are not in the paper, though everything that is reported in Table  
**    1 is computed 
*/

version 14
clear all
/*This is done once and for all including setting any parms used and controlling looping*/
set matsize 2000
set seed 12947
set more off
capture log close
log using bias.log, replace
/*the stuff below may change for different sets of runs*/
   matrix times = (100\75\50\30\20\10\7\5\3)
local ntimes=rowsof(times)
matrix groups = (20\50\100)
local ngroups=rowsof(groups)
matrix conterm = (0)
matrix corrterm =(2)
local ncorr = rowsof(corrterm)
  local nnovar=0 /*number of units where all y are 0*/
local  sims=	1000	/* number of simulations */
 matrix xco=(1)
local nx = rowsof(xco)
local ncons=rowsof(conterm)
local rowres=`ncons'*`ntimes'*`ngroups'*`ncorr'*`nx'
matrix define resnew=J(`rowres',7,0)

matrix colnames resnew = N T RMSLOGIT BIASLOGIT RMSCLOGIT BIASCLOGIT CORRTERM
/* 3 and 4 are RMS errors, 5 and 6 are signed errors, 7 is number of groups dropped */
noisily display "X is sum of SN mutliplied by xfecorr and the N fixed effects"
noisily display    " SIMS " `sims'
/*nothing changed below*/
local counter=0
  /* Now we loop over constants and T */

  forvalues xfn = 1/`nx' {
  forvalues consum=1/`ncorr' {
  forvalues groupnum = 1/`ngroups' {
forvalues connum=1/`ncons' {
forvalues  timenum = 1/`ntimes'  {
quietly {
local counter=`counter'+1
drop _all
set obs 10000
local xcoeff = xco[`xfn',1]
local xfecorr=corrterm[`consum',1]
local  time=times[`timenum',1]
 local n=groups[`groupnum',1]
local concoef= conterm[`connum',1] /* coefficient on constant term */
noisily display "sim number  " `counter' " n  " `n'  " t " `time'
local t=`time'			/* number of time periods */
local nomit=`nnovar'*`t'

local nt=`n'*`t'		
local n1=`n'+1
local n2=`n'+2
local t1=`t'+1


/*now we are going to generate the period and group variable*/
  
gen group = ceil(_n/`t') if _n<=`nt'
gen time = mod(_n-1,`t')+1 if _n<=`nt'
gen rangep = . /*80% range of probs*/
gen avt = . /*storing t ratios of clogit */
gen meanskew=. /*skewness in x*/
gen meankurt=. /*kurtosis in x*/
gen meanrsq = . /*rsq of x on fe1*/
gen errlogit = . /*absolute error in b for X coef of logit*/
    gen errlogittrue = . /*absolute error in b for X coef of logit w true dummy*/
gen errclogit = . /* absolute error in b for x coef of clogit*/
    gen coeflogit = . /*b for X coef of logit*/
        gen coefclogit = . /*  b for x coef of clogit*/
            gen errlogitnofe=.
gen coeflogitnofe=.
    gen meanp=. /*mean of the probs*/
      gen groupdrop = . /* groups dropped by clogit */
          gen perrlogit=.
gen perrlogitc=.
local i=1
while `i'<=`sims' {
/* now we are going to draw the true effects which do not vary by group */
generate fe=rnormal() if _n<=`nt'
generate fe1=fe
replace fe1 = fe1[_n-1] if group == group[_n-1]
/*we are going to give x some kurtosis by using a mild lognormal */
generate xstart= rnormal() if _n<=`nt'
 /*this is to get x corr with fe's* and then standardized */
egen x=std(`xfecorr'*xstart +fe1) if _n<=`nt'


/*this is to get the R^2 of x on fe1 */
  regress x fe1 in 1/`nt'
replace meanrsq=e(r2) in `i'
/*uncomment below if want to drop some groups with all y zero */
/*sort fefe group
local nomit=`n'*`nnovar'
replace fefe=-5 if _n<=`nomit' */

	/* draws Y from a Bernoulli distribution and estimates parameters */
gen p=	1/(1+exp(-(`concoef'+fe1+`xcoeff'*x))) if _n<= `nt'
                                            summarize p, detail
                       replace meanp = r(mean) in `i'
                       replace rangep = r(p90)-r(p10) in `i'
                     /*here is where we pick up the old sim code */
                       
	gen temp=uniform() if _n<=`nt'
	gen y=0
replace y=1 if temp<p



                    
capture clogit y x in 1/`nt', group(group) iterate(10)
            pause           
                        if _rc!=0 | e(ic)>9 {
                            noisily  display "bad clogit"
                            pause
                       capture drop temp   y   p   x fe fe1 xstart
                         continue
                       }
        replace errclogit=(_b[x]-`xcoeff')^2 in `i'
replace coefclogit=_b[x] in `i'
   local tc=_b[x]
                  replace groupdrop = e(N_group_drop) in `i'
        replace groupdrop = 0 if groupdrop == . in `i'

capture 	logit y i.group  x in 1/`nt', iterate(10)
       
if _rc!=0 | e(ic)>9  {
                         noisily display "bad logit"
                         capture drop temp  y   p   x fe fe1 xstart
                         continue
                       }
        replace errlogit=(_b[x]-`xcoeff')^2 in `i'
replace coeflogit=_b[x] in `i'


/*
constraint 1 x=`tc'

logit y i.group x in 1/`nt', constraint(1)
                       
                     predict logprobc
gen mslogprobc=(p-logprobc)^2
summar mslogprobc
replace perrlogitc=sqrt(r(mean)) in `i'

capture drop logprob logprobc mslogprob mslogprobc 

/*capture 	logit y fe1 x in 1/`nt', nocons
replace errlogittrue=(_b[x]-`xcoeff')^2 in `i'*/
*/
capture drop temp  y    p  x fe fe1 xstart

	local i=`i'+1

      }
           
/*this is the end of a single sim */



/*done after all the sims for a given T and concoef are done */

matrix resnew[`counter',1]=`n'  /*number of groups*/
  matrix resnew[`counter',7]=meanrsq[1] /*corrterm*/
  

matrix resnew[`counter',2]=`time' /*number of time periods*/

 





summarize errlogit
matrix resnew[`counter',3]=sqrt(r(mean))
summarize coeflogit
matrix resnew[`counter', 4] = r(mean)-`xcoeff'
summarize errclogit
matrix resnew[`counter',5]=sqrt(r(mean))
summarize coefclogit
matrix resnew[`counter', 6] = r(mean)-`xcoeff'

noisily {
 
matrix list  resnew
}

}  /*quietly*/
}  /*xcoef*/
} /* corr*/
} /* groups */
} /* constant */
} /* time */
/*Done once at the end after all looping is done*/
noisily {
 
matrix list  resnew
}

capture erase logitrms2.csv
mat2csv, matrix(resnew) saving(logitrms2)

